% Test of detection basing on reversed image after background subtration  %
% Example
% >> data = pathology_init_data(name);
% >> fungi_fluocell_detect(data);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%Lexie and Kathy Lu since 1/29/2015%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% input the information of the image, i.e. the path, name, and format

function [num_objects, fg_prop, colorvec, proj_vectors, cropped_images, orientations, boundary_roughness_vector] = fungi_fluocell_detect(data)
%    p = 'C:/Users/wang-lab/Desktop/pathology/data/0622/fungi/';
    p = data.path;
    file_name = data.file;
    fungi_intensity_th = data.intensity_th;
    % minimal and maximal area for fungi in image
    min_area = data.min_area; 
    max_area = 20 * min_area;

    %% load images and preproces the images
    im = imread([p, file_name]);
    figure; imagesc(im);
    % reverse the image before background subtraction
    rev_im = 255 * uint8(ones(size(im))) - im;
    % subtract background on counter staining
    [~, name, ~] = fileparts(file_name);
    file_bg = strcat(p, 'output/background_', name, '.mat');
    file_bg_value = strcat(p, 'output/background_value.mat');
    if ~strcmp(data.type, 'Scanscope') || ~exist(file_bg_value, 'file'),
        [bw, poly] = get_background(im, file_bg);
        figure; imshow(rev_im);
        hold on; plot(poly(:,1), poly(:,2), 'b', 'LineWidth',1.5);
        bg_mean = zeros(3, 1);
        for i = 1:3,
            bg_mean(i) = compute_average_value(rev_im(:,:,i) ,bw);
        end;
        if strcmp(data.type, 'Scanscope'),
            save(file_bg_value, 'bg_mean');
        end;
    else % data.type is Scancope, and exist file_bg_value
        temp = load(file_bg_value); 
        bg_mean = temp.bg_mean;
    end;
    % Background subtraction followed by converting RGB to Index
    for i = 1:3,
        rev_im(:, :, i) = uint8(rev_im(:, :, i) - bg_mean(i));
    end;
    % im_prep = sum(uint16(rev_im),3); gives the same result
    % since the sum function converts the variables to double.
    % The double type of rev_im also works better than uint16 type in
    % the my_graythresh function
    im_prep = sum(double(rev_im), 3);

    %% Detect the fungi/objects
    % find the global threshold of the image 
    % The maximum of an automatic threshold and an intensity threshold is used. 
    % Objects less than min_area and larger than max_area are both
    % removed.
    level = my_graythresh(im_prep);
    fg_bw = im2bw(im_prep, level) & (im_prep > fungi_intensity_th);
    % define a minimal area for fungi and get rid of staining dots
    temp = bwareaopen(fg_bw, min_area); clear fg_bw;
    fg_bw = temp;
    cc = bwconncomp(fg_bw, 8);
    fg_prop = regionprops(cc,'Area', 'PixelIdxList');
    fg_keep = (cat(1, fg_prop.Area) < max_area);
    pixel_index_keep = cat(1, fg_prop(fg_keep).PixelIdxList);
    temp = false(size(fg_bw));
    temp(pixel_index_keep) = true; clear fg_bw;
    fg_bw = temp; clear temp;
    clear temp cc fg_prop fg_keep clear pixel_index_keep;
    
    %% Calculate the characteristics of the detected fungi,
    % such as area, color, aspect ratio, etc.
    % The color is the color of the original image.
    
    % Use the watershed method to separate the bulk detection
    % The average color is the original color. 
    temp = my_WatershedSeg(fg_bw, data.type); clear fg_bw;
    fg_bw = temp; clear temp;
    % get the boundary of the detections
    % Find the accurate detection and get the mean value of each channel
    cc = bwconncomp(fg_bw, 8);
    num_objects = cc.NumObjects;
    avg_color = zeros(num_objects, 3); % for original 3-d image
    %avg_color = zeros(num_objects, 1); % for reversed background subtaction image
%    fg_mask = cell(num_objects, 1);
    for i = 1 : num_objects,
%         % Kathy : This is not efficient, I made the modification below
%         % --- 03/24/2015
%         % It is more efficient to use PixelIdexList to index into the image
%         % matrix
%         im_plain = false(size(im_prep));
%         im_plain(cc.PixelIdxList{i}) = true;
%         fg_mask{i} = uint8(im_plain); % VERY INEFFICIENT
%         for j = 1:size(avg_color,2),
%             %get the mean value on each channel
%             avg_color(i, j) = sum(sum(fg_mask{i} .* im(:,:,j))) / length(cc.PixelIdxList{i}); %based on original image           
%         end
%
        pixel_index = cc.PixelIdxList{i};
        %get the mean value on each channel
        for j = 1:size(avg_color, 2),
            temp = im(:, : , j);
            avg_color(i, j) = sum(temp(pixel_index))/length(pixel_index); 
            clear temp;
        end;
        clear pixel_index;
    end; 
    
    % Get the properties of each detection
    fg_prop = regionprops(cc, im_prep, 'PixelIdxList','Eccentricity','Solidity','Area');
    % Calculate and display the eccentricity of the detected objects
%     fg_eccentricity = zeros(size(fg_bw));
%     for i = 1:num_objects,
%         fg_eccentricity(fg_prop(i).PixelIdxList) = fg_prop(i).Eccentricity;
%     end;
%     figure; imagesc(fg_eccentricity);
%
%     % Solidity Area/ConvectArea
%     fg_solidity = zeros(size(fg_bw));
%     for i = 1:num_objects,
%         fg_solidity(cc.PixelIdxList{i}) = fg_prop(i).Solidity;
%     end;
%    figure; imagesc(fg_solidity);
    % Boundary roughness
    fg_bd = bwboundaries(fg_bw, 'noholes');
    fg_bd_roughness = zeros(size(fg_bw));
    boundary_roughness_vector = zeros(num_objects, 1);
    for i = 1:num_objects,
        % first direvative along the x and y, since fb_bd{i} is nx2
        boundary_roughness_vector(i) = sum(sum(abs(diff(fg_bd{i}, 1, 1))))/...
            sqrt(fg_prop(i).Area);
        fg_bd_roughness(cc.PixelIdxList{i}) = boundary_roughness_vector(i);
    end;
    figure; imagesc(fg_bd_roughness); title('Boundary Roughness');
    
    % the following assumes that there are actually detections
    if(not(isempty(fg_prop)))
        fg_keep = (cat(1, fg_prop.Eccentricity) <= data.eccentricity_th) & ...
            (cat(1, fg_prop.Solidity) >= data.solidity_th)& ... 
            boundary_roughness_vector <= data.roughness_th;
        pixel_index_keep = cat(1, fg_prop(fg_keep).PixelIdxList);
        temp = false(size(fg_bw));
        temp(pixel_index_keep) = true; clear fg_bw;
        fg_bw = temp; clear temp;
        figure; imagesc(fg_bw);
        clear fg_prop fg_bd num_objects;
        fg_prop = regionprops(fg_bw, im_prep, 'all');
        fg_bd = bwboundaries(fg_bw, 'noholes');
        num_objects = length(fg_prop);
    else % in the case that there are no detections
       clear fg_prop fg_bd num_objects;
       fg_prop = regionprops(0, 'all'); % we still want this to have the regionprops struct
       num_objects = 0;
       fg_bd = bwboundaries(fg_bw, 'noholes');
    end
    
    % Projection feature extraction
%     fg_bwlabel = bwlabel(fg_bw);
%     figure; imagesc(fg_bwlabel);
%     orientationv = -[fg_prop.Orientation]';
%     num_fungi = length(orientationv);
    cropped_images = cell(num_objects, 1);
    proj_vectors = cell(num_objects, 1);
    %%% This part broke and was commented Kathy - 03/11/2015
%     for i = 1:num_fungi
%         cropped_images{i} = cropper(im, fg_bwlabel, fg_bd{i}, i);
%     end;
%     for i = 1:num_fungi
%         proj_vectors{i} = sum(imrotate(cropped_images{i}, orientationv(i)));
%     end;
    %% plot a graph of average color of each item based 
%     figure;
%     scatter3(avg_color_list(:,1), avg_color_list(:,2), avg_color_list(:,3), 'fill');
%     hold on
%     scatter3(avg_color(:,1), avg_color(:,2), avg_color(:,3), 'fill', 'r');
%     legend('Jonathan', 'Lexie');
% plot the average color based on the reversed background subtraction image
    % figure; plot(1: length(avg_color), avg_color, 'b*');
    figure; plot(1: length(avg_color), avg_color(:,1), 'b*');
    xlabel('Index'); ylabel('Average Color (Red Channel)');
    %% show the final detection
% image on colour
    figure; imagesc(im); hold on; title('Final Detection');
    orange = [216, 41, 0]/255;
    
    for i = 1:num_objects,
        plot(fg_bd{i}(:,2), fg_bd{i}(:,1), 'color', orange, 'LineWidth', 1.5);
    end;
% % label each item marked as postive, respective to the index of avg_color 
%     for i = 1:size(avg_color,1),
%         s = sprintf('%d', i);
%         text(floor(fg_prop(i).Centroid(1)), floor(fg_prop(i).Centroid(2)), s, 'BackgroundColor',[.7 .9 .7]);
%     end;
% image on grayscale
    figure; imagesc(im_prep); hold on;
    colormap(gray);

    for i = 1:num_objects,
        plot(fg_bd{i}(:,2), fg_bd{i}(:,1), 'color', orange, 'LineWidth', 1.5);
    end;
  colorvec = avg_color;
  orientations = regionprops(fg_bw, 'Orientation');
return;
